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We study the phase structure of QCD at high temperature and density by lattice QCD simulations 
adopting a histogram method. We try to solve the problems which arise in the numerical study 
of the finite density QCD, focusing on the probabiUty distribution function (histogram). As a 
first step, we investigate the quark mass dependence and the chemical potential dependence of 
the probability distribution function as a function of the Polyakov loop when all quark masses 
are sufficiently large, and study the properties of the distribution function. The effect from the 
complex phase of the quark determinant is estimated explicitly. The shape of the distribution 
function changes with the quark mass and the chemical potential. Through the shape of the 
distribution, the critical surface which separates the first order transition and crossover regions in 
the heavy quark region is determined for the 2+1-flavor case. 
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1. Histogram method 

Not only the temperature (T) and chemical potential (/i) but also the quark masses are im- 
portant to understand the properties of QCD phase transition. In fact, the structure of the phase 
boundary in T — /i phase diagram is rather sensitive to the value of the strange quark mass. Recent 
lattice QCD simulations suggest that, for physical quai^k masses, finite T transition is crossover at 
zero jj., while it becomes first order for sufficiently large jj.. Identifying the critical point separating 
crossover and first order is one of the most challenging topics in lattice QCD simulations and in 
heavy-ion experiments. Probability distribution function or the histogram of the order parameter 
provides us with an important clue to identify such point in numerical simulations: In the case of 
the first order transition, different phases coexist at the transition point, so that the probability dis- 
tribution function has multiple peaks. On the other hand, in the case of crossover, such phenomena 
does not take place. Therefore, the nature of the transition can be identified through the shape of 
the distribution function. 

In this paper, we study the boundary of the first order transition region in QCD in the case 
when quarks are all heavy. We determine the boundary as function of the chemical potential jJ. 
by measuring histograms. Although this boundary in the heavy quark region is irrelevant to the 
boundary near the physical point, this provides us with a good testing and developing ground for 
the method, because the computational burden is much lighter. 

Selecting a physical quantity X, we calculate the probability distribution function defined by 



r Nf 
w{X,p,Kf,^f) = / m&Y^V S{X-X) g-^^-^s = 5{X-X) e-^^ ]^detM(K:j,Ai/) 

J J /=i 



\^ = l ^ ' / (X fixed;i3) 



where Sg, Sq and detM are the gauge action, the quark action and the quark determinant, respec- 
tively. Kf is the hopping parameter for the flavor quark mass. j8 = is the gauge coupling, 
and A^f is the number of flavors. (•• fixed;/?) = {' " ~ ^)) p / {^{^ ~ ^)) p means the expec- 
tation value measured with fixing the operator X at j8 in quenched simulations, Kf = /ly = 0. The 
expectation value in the right hand side is the ratio of w(X,j3, K.f,}JLf) and w(X,j8,0,0). However, 
the calculation of detM is usually difficult. We perform the hopping parameter expansion and 
compute the quark determinant in the leading order of the expansion |jl]]. 



detM(K-,jU) _ 
detM(0,0) 



288AfsiteK-'^^ + SAfs'l^'+^K^' jcosh (^^) + /sinh ^ij 



(1.2) 



for the standard Wilson quark action, where detM(0,0) = 1. The number of sites is A'^site = A'^^^ x 
Nf The quark determinant is simply given by the average plaquette operator P and the real and 
imaginary parts of the Polyakov loop operator, Q. = + iQ.\. Because the critical K is very small, 
at least for Nt = 4, this approximation can be justified for the determination of the critical K. 

In this calculation, it is essential to perform simulations at several simulation points and to 
combine these data by the multi -point reweighting method [^. Since values of the most observables 
distribute in a narrow range during one Monte-Carlo simulation, it is difficult to investigate the 
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shape of the distribution in a wide range. We thus combine several simulations. The expectation 
value of a operator X at jS is computed by simulations at j8; with the number of configuration A^, 
for / = 1, • • • ,Nsp using the following equation for the case of the plaquette action and degenerate 
A'^f-flavor, 

where the weight factor G{P) is 

G{P) = z , (1.4) 

and (• • •)aii means the average over all configurations generated at all j8; with K = fj. =0. The par- 
tition functions iF(j8,) ai^e pai^ameters in this method and are determined by solving a consistency 
condition: i2°(j3,) ^j^jj ^^^^ ^ G{P) , numerically for each / = 1, • • • ,A'^sp> except for an overall 
normalization constant. ^^^f j means the sum of configurations at all /3;. (See the appendix A 
in Ref. ^ for details.) Note that this method enables us to change K and /3 continuously. 



2. Polyakov loop distribution function at zero density 

The most important observable near the transition point in the heavy quark region is the 
Polyakov loop, which is the order parameter of the deconfinement transition. We analyze the data 
obtained at 5 simulation points, j8 = 5.68 - 5.70, in the quenched simulations with the plaquette 
gauge action on a 24-^ x 4 lattice [|l|]. Figure [l] is the result of the Polyakov loop susceptibility, 
Xq = N^{{D. - as a function of K^' and the effective j3 defined as /3* = j3 +4SNfK^ with 



A^, = 4, computed at = using Eq. (1.3). Because the plaquette action is Sg = —SNsiteliP, the 
plaquette term in Eq. (1.2) can be absorbed into the gauge action by defining /3* and the analysis be- 
comes simpler. Owing to the multi-point reweighting method, the susceptibility can be calculated 
in a wide range of jS and fc. We define the transition point as the peak position of Xa.- 
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Figure 3: fc-dependence of the probability distribution of the Polyakov loop at the transition point in the 
complex plane for A^f = 2. The value of k'^ is shown in the upper right of the figure. 

We then measure the distribution function of the absolute value of the Polyakov loop at the 
transition point for the case of 2-flavor QCD at jU = 0. Expectation values with fixing the value 
of the Polyakov loop are computed using the delta function approximated by a Gaussian function, 
5{x) exp[—{x/A)^]/{A^/7i), where A = 0.005 is adopted consulting the resolution and the sta- 
tistical error. We plot the effective potential, Veff(|^^|) = — lnw(|n|), for several values of in 
Fig. ^. j8 is adjusted to the peak position of Xo. at each K. The value of Veff(|n|) is normalized at 
|n| = 0.01. This figure shows that the shape of Veff(|^^|) is double-well type at = 0, indicating 
the first order transition, and the shape changes gradually as increasing K. It becomes single-well 
around ~ 0.00002, suggesting the first order transition changes to crossover. The critical value 
of K has been determined by measuring the distribution function of the average plaquette in Ref . [|l|] 
with the same configurations. The result is K^p = 0.0658(3) (^j^ ) for N{ = 2. Hence, the results of 
Kcp from the plaquette and Polyakov loop effective potentials are consistent with each other. 

We moreover calculate the distribution function of the complex Polyakov loop in the complex 
plane (Hrj^i) at the phase transition point, which is shown in Fig. |3|for 2-flavor QCD. The well- 
known Z(3) symmetric 4 peak structure is observed at k: = 0, and the 2 peaks in the negative 
region become smaller as increasing K. Then, the remaining 2 peaks are getting closer with K, and 
the distribution becomes a single peak around K = 0.00002. These figures illustrate how the Z(3) 
symmetric quenched QCD changes to full QCD. 

3. Complex phase and distribution function at finite density 

The histogram method is powerful in particular with the presence of the chemical potential 
/I. Direct simulations by the Monte Carlo method cannot be performed at finite chemical potential 
because the quark determinant is complex. An approach to simulate finite density QCD is to 
combine the reweighting method and simulations with the complex phase of the quark determinant 
suppressed, which are called phase-quenched simulations. The distribution function for the real 
part of Polyakov loop, Hr, is a good example to explain the contribution from the complex phase 
and the phase-quenched part in the distribution function at finite density. 
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Figure 4: The solid lines are yeff(nR) at = 
for each k^. VeffC^^R) at finite K''*cosh(;U/r) is 
between the solid line and the dashed line. 
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Figure 5: The average of the complex phase fac- 
tor and the 2"'', 4* and 6* order cumulants calcu- 
lated with fixed Q.r at K-'*sinh(^/r) w 0.00002. 



We calculate the distribution function in heavy quark QCD for the degenerate Nf = 2 standard 
Wilson case. Using the hopping parameter expansion, the distribution function can be factorized 
into the phase factor and the phase-quenched part. 



w(nR,/3,JC,At) = / W5(ilR-aR)/^^-'=/5^(detM(jc,At)) 



:w(aR,j8,0,0)(e 



,288A'si,eWfK-'P 



exp 



3N^2^'+'-NfK^' jcosh (^jj CIr + /sinh 



w 



(aR,/3*,0,0)exp [SN^l'^'+^NfK^' cosh Qr 



{£^R;i3*,0) 



(3.1) 



where 6 is the phase of the quark determinant: 

= 3N^2^'+'-NfK^' sinh(Ai/r) ili, 



(3.2) 



and the part in front of the phase average is the distribution function in the phase-quenched theory. 
The plaquette term is absorbed by Sg shifting j3 to /3* = j3 +48^^:^. {■ ■ •) {ag_;ii.K) means the 
expectation value at (j8, k) with fixing 

The phase-quenched part of w{D.r,[5, k:,}Jl) can be obtained from that at /i = simply by re- 
placing K"^' by K^' cosh(/i/r), since the distribution function at = is given by w{Q.r,^* ,Q,Q) 
X exp[3A'^^^2'^'+^K■^'^lR]. Therefore, the critical value in the phase-quenched theory is given 

by K^'(O) = K"^'(/x)cosh(/i/r). Moreover, adopting K"'^' cosh(/i/r) as basic parameter to in- 
vestigate the critical point, the magnitude of the phase is limited for each K"'^' cosh(jU/r) be- 
cause k:'^' sinh(/i/r) in d is always smaller than K:'^'cosh(/i/r). The effective potential of Hr, 
Veff(^R) = — InH'(nR), at the transition point for /i = is plotted in Fig. Qby the sohd lines for 
each K^< , and is equal to the phase-quenched distribution function for each k^' cosh (^u/r). This 
Veff is normalized at Hr = 0. 

Next, we calculate the phase factor, (^'^)(£1r;/3*,o)- If ^'^ changes its sign frequently, the statis- 
tical error becomes larger than the expectation value, causing the sign problem. To avoid the sign 
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problem, we evaluate the phase factor by the cumulant expansion |5|] : 

(3.3) 



= exp 



. {B\ i{d\ {e% i{e% {6% 

'^^^'-—1 3r+^+^^ 6r+' 



where (0"), is the order cumulant: {9^), = {d^)^ar,fi*fl), {e^)c = (0^)(nR;r,o) -3(0^>ffi^.p.,o)' 
{e')c = (0')(nK;r,o) - 15(0^)(nR;/3*,o)(0'){a,;r,o) +30(02)3^^.^^^^^, .... The key point of this 
method is that {d")c = for any odd n due to the symmetry under 6^—6, and thus the complex 
phase can be omitted from this equation. This implies that {e'^) is guaranteed to be real and positive 
and the sign problem is resolved once the cumulant expansion converges. Another important point 
is that 6 is given by the average of the Polyakov loop. When the correlation length is finite, the 
phase can be written as a summation of local contributions G = dx with almost independent dx- 
The phase average is then 



n^"' =-pII^(^">j- (3-4) 



This suggests that all cumulants (0")c ~ 'Lx{^x)c increase in proportion to the volume as the vol- 
ume increases. Only for such a case, the effective potential Veff can be well-defined though the 
phase-quenched effective potential Vq with Veff = Vq — ln{e'^) = Vq — L«'"(^")c-/"' in large 
volume limit, since Veff and Vq are both in proportion to the volume. 

We plot the results of {d")c/n\ in Fig. | for fc^' sinh(;U/r) = 0.00002 and j3* = 5.69. The 
black, blue and green lines are the results for « = 2,4 and 6, respectively. The fourth and sixth 
order cumulants are very small in comparison to the second order for this k. The red line is 
— In(e'^)(j2j^.p* 0)' which is almost indistinguishable from the second order cumulant. The contribu- 
tion from the fourth and sixth orders becomes visible at small as K"'^' sinh(/x/r) increases. How- 
ever, for the determination of the critical point, cosh(/i/r) 0.00002 in the phase-quenched 
theory, the region at K^' sinh(jU/r) < 0.00002 is important because cosh(/x/r) > sinh(/i/r). This 
figure thus indicates that the phase average is well-approximated by the second order cumulant 
around the critical K. The results of the effective potentials including the effect from the phase 
factor are shown by the dashed lines in Fig. ^ for each k'^' cosh(jLt/r). In this figure, the phase 
factor is estimated by the second order cumulant at /i/r = oo, i.e. sinh(/i/r) = cosh(jLt/r). The 
results at finite jj. are between the solid and dashed lines. We find that the contribution from the 
phase, —ln{e'^), is quite small except at small and the phase factor does not affect Veff in the 
region of relevant for the determination of the critical point. This means that the contribution 
from the complex phase to the location of the critical point is quite small on our 24^ x 4 lattice. 

Neglecting the effect of the phase factor, it is easy to determine the critical point for the A^f = 
2+1 case because the difference from the Nf = 2 case is just to replace Ik'^' by 2k^^ + K^'. We 
thus find that the critical (K"ud, K"s) is given by 

2Kf^(At)cosh(Atud/r) + Kf'(At)cosh(Ats/r) =2[K^f=2(o)]iv,^ (35) 

where j4'=^(0) = 0.0658(3)(lii ) for Nt = 4 [0]. The critical lines in the K plane for up, down and 
strange are drown in Fig. ^ for the cases of Hud/T = jJ-a/T = 0-10 (left) and jJ-ud/T = 0-10 and 
Ats/r = (right). 
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Figure 6: Critical surface separating the first order transition and crossover regions in the heavy-quark 
region. Left: The case fi,, — jJ-d ~ fix = Right: The case that may be reahzed in heavy ion colUsions; 
liu = Hd = ^„d and pLs = 0. 

4. Summary 

We have studied the phase structure of QCD at nonzero chemical potential pL in the heavy quark 
region, highlighting the properties of the probability distribution function of the Polyakov loop Q.. 
The shape of the effective potential defined by the distribution function changes with the quark 
mass and the chemical potential. The multi-point reweighting technique enables us to obtain the 
distribution function in a wide range of the coupling parameters. We have shown that the effective 
potential provides us with an intuitive and powerful way to investigate the fate of first order phase 
transitions. Through the shape of the potential, the critical surface where the first order deconfining 
transition in the heavy quark limit terminates is determined for the 2-1-1 -flavor case. The effect from 
the complex phase of the quark determinant has been estimated explicitly, and is found to be quite 
small around the critical point for any chemical potential in the heavy quark region. On the other 
hand, the effect from the complex phase must be important in the light quark region. An attempt to 
study finite density QCD at light quark masses by combining phase-quenched simulations and the 
reweighting technique is reported in 
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